
use book-basedata-replication, clear

quietly stcox med2 med2_t prevcris2 viol2 crisdur2 jointdem victory2 contig2 /* 
*/		prevcris2_t viol2_t crisdur2_t jointdem_t victory2_t contig2_t /* 
*/		if _t<3650,  strata(order5) cluster(dyadno) efron nohr


*calculate mean of each variable
foreach var of varlist med2 prevcris2 viol2 crisdur2 jointdem victory2 contig2 {
quietly sum `var' if e(sample)
local s_`var'=r(mean)
display `s_`var''
}


scurve_tvc if _t<3650, gen(base_S_mean) at(med2 0 prevcris2 `s_prevcris2' /* 
*/		viol2 `s_viol2' crisdur2 `s_crisdur2' jointdem `s_jointdem' /* 
*/		victory2 `s_victory2' contig2 `s_contig2') tvc(med2) texp(_t) /*
*/ 		strata(order5) ties(efron)


scurve_tvc if _t<3650, gen(med_S_mean) at(med2 1 prevcris2 `s_prevcris2' /* 
*/		viol2 `s_viol2' crisdur2 `s_crisdur2' jointdem `s_jointdem' /* 
*/		victory2 `s_victory2' contig2 `s_contig2') tvc(med2) texp(_t) /*
*/ 		strata(order5) ties(efron)




*label variables
label var base_S_mean1 "No mediation"
label var med_S_mean1 "Mediation"



*estimate Royston-Parmar model for the same model specification
xi: stpm2 med2 prevcris2 viol2 crisdur2 jointdem victory2 contig2 i.order5 /* 
*/		if _t<3650, scale(hazard) df(5) tvc(med2 i.order5) dftvc(med2:1 5)

*predict survival function at mean covariate values for mediated/unmediated case
predict s0_rp, survival at(med2 0 prevcris2 `s_prevcris2' viol2 `s_viol2' /* 
*/		crisdur2 `s_crisdur2' jointdem `s_jointdem' victory2 `s_victory2' /* 
*/		contig2 `s_contig2') time(_t) ci
predict s1_rp, survival at(med2 1 prevcris2 `s_prevcris2' viol2 `s_viol2' /* 
*/		crisdur2 `s_crisdur2' jointdem `s_jointdem' victory2 `s_victory2' /* 
*/		contig2 `s_contig2') time(_t) ci
predict sdiff_rp, sdiff2(med2 0 prevcris2 `s_prevcris2' viol2 `s_viol2' /* 
*/		crisdur2 `s_crisdur2' jointdem `s_jointdem' victory2 `s_victory2' /* 
*/		contig2 `s_contig2') sdiff1(med2 1 prevcris2 `s_prevcris2' /* 
*/		viol2 `s_viol2' crisdur2 `s_crisdur2' jointdem `s_jointdem' /* 
*/		victory2 `s_victory2' contig2 `s_contig2') time(_t) ci



*label variables
label var s0_rp "No mediation"
label var s1_rp "Mediation"


*graph results for each stratum
twoway line base_S_mean1 med_S_mean1 _tscurve if _tscurve<3650, sort c(J J) /* 
*/		lcol(gs7 gs10) || /* 
*/		line s0_rp s1_rp _t if _t<3650 & order==1, sort lpattern(shortdash shortdash) /* 
*/		legend(order(1 2) size(small) region(lwidth(none))) lcol(gs1 gs6) /* 
*/		saving(order1, replace) t2title("1 previous crisis") ylab(0(.2)1)
forvalues t=2/4 {
twoway line base_S_mean`t' med_S_mean`t' _tscurve if _tscurve<3650, sort c(J J) /* 
*/		lcol(gs7 gs10) || /* 
*/		line s0_rp s1_rp _t if _t<3650 & order==`t', sort lpattern(shortdash shortdash) /* 
*/		lcol(gs1 gs6) saving(order`t', replace) legend(off)  /* 
*/		t2title("`t' previous crises") ylab(0(.2)1)
}
twoway line base_S_mean5 med_S_mean5 _tscurve if _tscurve<3650, sort c(J J) /* 
*/		lcol(gs7 gs10) || /* 
*/		line s0_rp s1_rp _t if _t<3650 & order==5, sort lpattern(shortdash shortdash)  /* 
*/		lcol(gs1 gs6) saving(order5, replace) legend(off)  /* 
*/		t2title("5 or more previous crises") ylab(0(.2)1)
